Feature/local background subtraction - #131

Hi Manuel,

thanks for adding the feature. For traceability can you maybe elaborate the feature and model behind that here?

Thank you!

Hi Kristian,

sure — here is the feature and the model behind it.

What the feature adds. Step 3 computes one additional quantity per LED and image, bgsub_sum_col_value: the local background $\hat B$ is estimated as the median of the positive pixels in the outer 2-px border ring of the LED search area, and then subtracted from every pixel before summing:

$$\hat B = \mathrm{median}\lbrace p_i \in \mathrm{ring} \mid p_i > 0 \rbrace, \qquad S_{\mathrm{bgsub}} = \sum_{\mathrm{box}} \max(p_i - \hat B,\ 0)$$

Median → robust against LED tails and hot pixels; positive-only → raw Bayer arrays mask foreign-channel pixels with exactly 0, which would otherwise pull the median to ~0; clip at 0 → masked pixels stay untouched. The analysis uses it only when reference_property = bgsub_sum_col_val is set in config_analysis.ini — the default pipeline is unchanged (verified bit-identical by the new acceptance suite 07, which also covers reading pre-existing 5-column step-3 files).

The model behind it. The search-area sum contains, besides the LED, an additive background (ambient/stray light plus a small estimator noise floor) that does not dim with smoke. Because this floor $F$ is present in the reference image and the measurement image, the measured transmission becomes (exactly, for a constant floor):

$$T_{\mathrm{meas}} = \frac{S_0,T + F}{S_0 + F} = T + f,(1-T), \qquad f = \frac{F}{S_0+F}$$

which propagates to the extinction coefficient as $\Delta k/k \approx -f,(e^{\tau}-1)/\tau$: the bias is always negative (extinction underestimated) and grows with optical depth. The ring measures $\hat B$ locally and per image, so the correction itself does not rely on the floor being constant in time or on any assumption about its origin.

Empirical basis (Hasenheide data, camera F1L). The LED-free border ring falls only to 0.86 of its reference value while the LED signal falls to 0.51 — direct evidence of light in the box that the smoke does not dim. The floor fraction measured independently from the ring time series ($f = 9.4\ %$ for a dim probe channel) predicts a $-11\ %$ bias via the formula above; $-12.6\ %$ was observed. For the bright production channels $f \approx 1.4\ %$, i.e. a $-2$ to $-4\ %$ effect. In a synthetic-twin test with a constant floor injected into the acceptance-test images, sum_col_val follows the formula prediction within 1 pp across $f = 0\ldots 32\ %$ while bgsub_sum_col_val stays within $+0.2\ %$ of the known ground truth.

Hi Kristian,

two test-only commits were added, so the step 3 output and bgsub_sum_col_val itself are unchanged:

• 7b65569: the bit-identity test in suite 07 never read the stripped step 3 files. read_hdf only builds all_parameters.h5 if it does not exist, and the binary from the previous suites was reused. The test now deletes the binary before each analysis run.
• 630fcec: the acceptance test data are generated with a fixed seed. The rmse of bgsub on the clean images varied between runs, from 0.039 to 0.066 against the bound of 0.06.

An update on the model, from the Hasenheide 2.0 data:

• With the small search areas (16–20 px), the 2-px border ring lies in the halo of the LED itself. The ring median was 163 / 358 counts, compared with 6 / 10 counts 25–60 px next to the strip. The halo scales with the LED and largely cancels in I/I0. So bgsub removes a mixture of halo and additive pedestal rather than the pedestal alone.
• In practice this works out well in weakly sooting tests. There, bgsub and a pedestal estimate from side bands next to the strip agree within 0–2 percentage points in LEDSA/LEM. Under dense smoke they differ: LEDSA/LEM on a plateau was 0.80–0.85 raw, 0.90–0.94 with bgsub and 0.96–1.03 with the side-band estimate.
• The "median of the pixels > 0" also overestimates a pedestal of a few counts. Many channel pixels are noise clipped at 0 after the black level subtraction, 35–45 % under smoke.

In my fork I built opt-in follow-ups on top of this PR. The defaults stay bit-identical:

• a side-band pedestal (sidebg_sum_col_val),
• a ring variant with a mean over the channel pixels including zeros (bgsub2_sum_col_val),
• a per-ray validity marking (-rv),
• an automatic reference image selection,
• a few small fixes.

If you are interested, I would offer them as separate PRs once #130 and #131 are merged. Would you merge these two as they are now?

Best,
Manuel

Hi Manuel,

thanks for the detailed write-up — the bias formula, the ring-vs-LED evidence from the Hasenheide data and the synthetic-twin check make the case well, and keeping it opt-in behind reference_property is the right call. The implementation of the estimator itself looks good to me: median over positive ring pixels, clip at 0, and the docstring explains why both choices are needed for raw Bayer arrays.

Two things before this can go in, one of them a blocker.

1. Blocker: -s3 followed by --analysis crashes

The fitted path is broken on this branch. Reproduced on 630fcec with the acceptance test data:

python -m ledsa -s1
python -m ledsa -s2
python -m ledsa -coord
python -m ledsa -s3          # fitted, not -s3_fast
python -m ledsa --analysis

ValueError: could not convert string '4.1210e+035.0521e+01' to float64 at row 0, column 6.
  File "ledsa/core/file_handling.py", line 96, in read_table
    data = np.loadtxt(filename, delimiter=delim, dtype=dtype)

The cause is not in your estimator, it is where the new column lands. LEDAnalysisData.get_data_string concatenates the main and the fit part without a separator:

out_str = self.get_main_data_string()
if self.fit_leds:
    out_str += self.get_fit_data_string()     # no ',' in between

get_main_data_string ends without a trailing comma and get_fit_data_string starts with led_center_x, so the last main field and led_center_x have always been glued together in fitted output. On master the last main field is max_col_value, an int, so the glued field reads e.g. 2325.0521e+01 — nonsense, but still a parseable float, and the resulting column count happens to match what _get_column_names returns. So master gets away with it: --analysis runs and the extinction coefficients are correct, because they are computed from sum_col_val, which sits before the damage. Everything from led_center_x onwards is silently shifted by one column, which is a separate bug I will open an issue for.

Your change moves bgsub_sum_col_value to the end of the main part, so the glued field becomes 4.1210e+035.0521e+01 — two floats in scientific notation with no delimiter, which numpy cannot parse. The latent bug turns into a hard failure.

Concretely, fitted step-3 output:

The fix is one comma in get_data_string, and then _get_column_names has to be corrected for the fitted format as well — see below.

2. _get_column_names and the (6, 20) check

if parameters.shape[1] in (6, 20):  # files written since the local background subtraction was added
    columns.append("bgsub_sum_col_val")

6 is right for -s3_fast. 20 never occurs: with the missing comma the fitted file has 19 columns, and once the comma is added it has 20 but the name list is still one short, because the existing list omits fit_time:

columns.extend(["x", "y", "dx", "dy", "A", "alpha", "wx", "wy", "fit_success", "fit_fun", "fit_nfev"])

That missing name is what makes the final parameters.shape[1] != len(columns) - 1 guard fail, which sends the fitted case to _get_old_columns. That returns [] for anything other than 15 or 4 columns, so _get_column_names ends up returning just ['width', 'height'].

I would rather not keep extending the positional heuristic — the next column added anywhere in the row breaks it again silently. Deriving the names from the exact column count reads better and fails loudly:

base_columns = ["img_id", "led_id", "led_array_id", "sum_col_val", "mean_col_val", "max_col_val"]
bgsub_columns = ["bgsub_sum_col_val"]
led_center_columns = ["led_center_x", "led_center_y"]
fit_columns = ["x", "y", "dx", "dy", "A", "alpha", "wx", "wy", "fit_success", "fit_fun", "fit_nfev", "fit_time"]

and then match parameters.shape[1] + 1 against the known combinations, falling back to _get_old_columns only for the genuinely old layouts and raising if nothing matches. Note 'img_id' is not in the file but prepended on load, which is where the - 1 in the current guard comes from; that off-by-one is easy to get wrong and worth a comment.

Your TestColumnNames covers the 6- and 5-column cases, which is why this stayed hidden — a fitted-format case there would have caught it. Could you add one?

I tried both changes together on your head: the fitted file then has 20 columns, -s3 followed by --analysis goes through, and your suites stay green (12/12 acceptance, 11/11 unit). Happy to send the patch if that saves you the round trip.

Thanks again — the feature itself I am happy with.

Hi Kristian,

thanks for the careful review and for tracking this down. You are right on both points. I reproduced the glued field on 630fcec: the fitted row has 19 columns and np.loadtxt fails on 4.1210e+045.0521e+01. Thanks also for the offer of the patch. I implemented it along the lines you sketched, so you don't need to send it:

• ba2cadb:
  • Delimiter: LEDAnalysisData.__str__ now puts a comma between the main and the fit part.
  • Column names: _get_column_names looks the names up by the exact number of columns, with BASE_COLUMNS, BGSUB_COLUMNS, LED_CENTER_COLUMNS and FIT_COLUMNS (now including fit_time) as you suggested. A comment explains the + 1 for img_id.
  • Known layouts: -s3_fast with 5 and 6 columns, and -s3 with 20 columns. The oldest layouts (4 and 15 columns) still go to _get_old_columns. Any other column count raises a ValueError naming the file and the count.
  • Tests: TestColumnNames now covers a fitted row written by LEDAnalysisData itself (it checks that max_col_val, led_center_x, A, fit_nfev and fit_time land in the right columns), the legacy fitted layout, the oldest layouts and unknown counts.
• 57bf5b0: a new acceptance test in suite 07 runs -s3 followed by --analysis with both reference properties. It fails on 630fcec (the binary cannot be built) and passes now. -s3 takes about two minutes on the test data.

One decision I'd like your view on: fitted files written by the current master (18 columns, with max_col_value and led_center_x glued). Raising for them would make existing fitted results unreadable, although sum_col_val and mean_col_val in them are correct. I kept reading this layout with the previous (shifted) names and print a warning that max_col_val and all fit columns are not valid and that step 3 should be run again for valid fit results. If you prefer a hard error there, it's a one-line change.

This also fixes the column shift of the fit columns that you wanted to open a separate issue for, so that issue may not be needed. If you open it anyway, I'll link it.

Unit tests: 18/18. Robot suites 02–07 pass, including the new test. S1 fails on Windows independently, see #135.

Best,
Manuel

Hi Manuel,

thanks for the quick turnaround — ba2cadb is what I had in mind, and the + 1 comment and the ValueError naming the file and the count are exactly the loud failure I was after. I re-ran your new test here: it reproduces the crash on 630fcec (could not convert string '4.1210e+035.0521e+01' to float64 at row 0, column 6) and passes at 57bf5b0, unit tests 18/18. atleast_2d=True taking the single-LED file with it is a nice extra.

On your question: keep the warning, please don't raise. sum_col_val and mean_col_val sit ahead of the glued field and are correct, and they are what the analysis actually uses, so the results from those files stay valid. A hard error would make all our pre-PR fitted data unreadable without buying any correctness.

In the same spirit, one addition: files written by this branch before the delimiter fix (19 columns, bgsub_sum_col_value glued to led_center_x) cannot be loaded at all — np.loadtxt fails before any column name is derived, and with silent=True the user sees only the bare ValueError. I ran into that while reproducing the bug, so everyone who tried the branch has such files. Could _get_column_names catch it and re-raise with the same hint as the 18-column warning, i.e. missing delimiter, re-run step 3?

master has moved since you branched

#137, #140 and #141 landed in the meantime, so please merge master in. Five things come with it:

• The conflict is only in ledsa/tests/AcceptanceTests/LedsaATestLibrary.py; file_handling.py, step_3_functions.py and ConfigDataAnalysis.py merge cleanly. We fixed the flaky noise twice, independently: master has a module-level NOISE_SEED = 20181127 used as np.random.default_rng(NOISE_SEED), while your 630fcec added a seed=1 argument with the legacy np.random.seed. Please keep master's version and drop the seed parameter — but your background= argument has no counterpart there and has to survive the merge (create_test_data, create_test_image, create_img_array).
• The thresholds in suite 07 need re-checking afterwards. default_rng(20181127) draws different noise than np.random.seed(1), so the synthetic ground truth changes, and rmse_sum < 0.05, rmse_bgsub < 0.06 and the rmse_sum > 0.05 assertion in the background test were all tuned on seed 1.
• Suite number. master now has 07_test_corrupted_images.robot, so your file should become 08_test_background_subtraction.robot.
• The encoding check in CI will fail on your new code. The Windows port added a step that greps ledsa/ for text I/O that does not name an encoding; master is clean, four new lines are not: LedsaATestLibrary.py:350,352 (open(path) and open(path, 'w') in strip_bgsub_column_from_led_position_files) and test_background_subtraction.py:102,133 (write_text(...)). encoding='utf-8' on all four.
• Windows is now gating (Let a failure of the Windows acceptance tests fail the workflow #140): the continue-on-error exemption is gone, since the Windows leg went green in make the data extraction and analysis run on Windows #139. The S1 failure you saw is most likely your branch predating that port, so it should go green with the merge — but anything left over now blocks instead of being tolerated.

One small note on the new test: on 630fcec it does not fail red. check_error_msg calls exit() as soon as the subprocess writes to stderr, so the whole robot run stops with "Execution stopped by user" instead of reporting one failed test. That is pre-existing and nothing for this PR — just so the commit message is not read as a clean red test by whoever looks next.

Once the merge is resolved and CI is green I am happy to approve. The separate issue for the fit column shift I will not open — ba2cadb covers it.

Best,
Kristian

Hi Kristian,

thanks, and for the clear guidance on the merge. All points are addressed:

• Merge (ecece77): the conflict in LedsaATestLibrary.py is resolved as you described.

  • master's NOISE_SEED with default_rng stays and the seed parameter is gone.
  • background= survives in create_test_data, create_test_image and create_img_array, as does the reference_property argument of load_extinction_coefficients_computed.
  • The suite is now 08_test_background_subtraction.robot.

• Encoding: encoding='utf-8' on the two open() calls in strip_bgsub_column_from_led_position_files (in the merge commit) and on the two write_text() calls in test_background_subtraction.py (573f57a). The grep of the CI step is clean on ledsa/.

• Unparsable 19-column files (6ae80d6): _get_column_names catches the ValueError of read_table and re-raises it with the file name and the hint that bgsub_sum_col_value and led_center_x are glued and step 3 has to be run again. A unit test covers it.

• Suite 08 needed a setup (ce7c23c): the teardown of the new suite 07 empties analysis/channel0, so the first two tests of the background suite found no step 3 files. The suite now runs -s3_fast in the work directory itself and removes a binary of an earlier run.

• Thresholds with NOISE_SEED, all within their bounds with margin:

  Test	rmse	Bound
  bgsub on clean images	0.045 / 0.034 / 0.029	< 0.06
  With background, sum_col_val	0.13 / 0.18 / 0.13	> 0.05
  With background, bgsub_sum_col_val	0.014 / 0.016 / 0.022	< 0.05
  Fitted output, sum_col_val	0.014–0.022	< 0.05
  Fitted output, bgsub_sum_col_val	0.029–0.045	< 0.06

• Your note on the new test is right, thanks: on 630fcec check_error_msg stops the whole run ("Execution stopped by user") instead of reporting one red test, so "fails" in the message of 57bf5b0 overstates it. I won't rewrite the pushed commit; this comment is the correction.

Checks before pushing:

• Windows, without UTF-8 mode (cp1252 locale): the Robot suite passes 17/17, including both S1 tests and the teardown. The unit tests pass 19/19.
• Default path, test data: -s3_fast followed by --analysis (linear and nonlinear, sum_col_val) gives bit-identical extinction coefficients on master and on this branch. The step 3 files are identical in their first 5 columns; the 6th is the new bgsub column.
• Default path, real data: the same on three raw images of our Hasenheide data, all three channels. There, bgsub_sum_col_val is also identical to the values we computed before your Windows port.
• Linux: I can only see the Linux leg in CI. Statically, the changed files parse with the Python 3.10 grammar, the committed blobs have LF line endings and git diff --check is clean.

One observation outside this PR, from checking the interplay with #137: _get_column_names always reads 1_led_positions.csv, and create_binary_data calls it before the loop that skips images without results. So if exactly image 1 could not be read, --analysis stops with FileNotFoundError: analysis/channel0/1_led_positions.csv not found instead of leaving the image out. This is the same on master and on this branch; I reproduced it by removing that file. Taking the first existing result file would fix it. I'm happy to do that in a separate PR if you like.

Best,
Manuel